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Abstract 


In this work laminar (that is, laminar liquid, laminar vapour) flow boiling of water as well 
as propane has been simulated qualitatively by a method known as the Coupled Map 
Lattice (CML) method. CML is a simple model introduced first in 1983 by Kunihiko 
Kaneko. The CML model has essential features of spatiotemporal chaos. Over the past 
two decades, studies of the CML have been expanding not only in the field of 
spatiotemporal chaos and pattern formation but also in the fields of biology, mathematics 
and engineering. In the area of pool boiling so far only four studies of CML have been 
reported. They are Yanagita (1992), Shoji and Tajima (1997), Shoji (1998) and 
Ghoshdastidar et al. (2004). So far no studies showing CML simulation of flow boiling 
have been formd in the open literature. The present work is a first attempt in this 
direction. 

A CML is a dynamical system with discrete-time, discrete-space and continuous 
states. It usually consists of dynamical elements on a lattice which interact (are ‘coupled’) 
with suitably chosen sets of other elements. The strategy of modelling dynamical 
phenomena in spatially extended systems by a CML model is based on the following 
steps. 

1 . Choose (set of) macroscopic variables on a lattice. 

2. Decompose the processes underlying the phenomena into independent components. 

3. Replace each component by a simple parallel dynamics on a lattice. 

4. Carry out each unit dynamics (or procedure) successively. 

The modelling by CML is based on the assumption that the flow boiling is governed 
by (a) nucleation from cavities on the heated surface (b) forced convection and (c) phase 
change in the fluid bulk and mixing. The macroscopic variable chosen is temperature. 
The problem is transient and marches to steady state solution. The forced convection part 
in CML part is solved using CFD techniques, assuming average properties of the fluid. 

The main observations from the results are the following, (i) The wall temperature 
does not affect the heat transfer coefficient and vapour fraction; (ii) The decrease in entry 
subcooling lowers the heat transfer coefficient near the tube entrance but does not 



influence the heat transfer coefficient in the latter part of the tube. Furthermore, vapor 
fraction is high near the tube entrance for lowAT^^^’s whereas the maximum vapour 

fraction is same for all ’s. (hi) The heat transfer coefficient is virtually independent 

of mass flow rate in the initial section of the tube where the nucleate boiling effect is 
dominant. On the other hand, in the latter part of the tube where the convection 
vapourization is significant, heat transfer coefficient falls with increase m mass flow rate. 
All these trends are similar for both water and propane. 

An exception for propane is noted where it is seen that fluid mean temperature is 
higher for = 20°C than for AT^^^ = 30°C in the section of the tube from z = 4.8m to 

6m. Another difference with respect to water is that the convective vapourization zone is 
much larger as compared to the zone where the nucleate boiling is dominant. In addition 
to the above, the drop in heat transfer coefficient at some location down the tube is not as 
marked as that for water. 
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CHAPTER 1 


INTRODUCTION 


1.1 Introductory Remarks 

In this work laminar (that is, laminar liquid, laminar vapour) flow boiling of water as 
well as propane has been simulated qualitatively by a method known as the Coupled 
Map Lattice (CML) method. CML is a simple model introduced first in 1983 by 
Kunihiko Kaneko. The CML model has essential features of spatiotemporal chaos. 
Over the past two decades, studies of the CML have been expanding not only in the 
field of spatiotemporal chaos and pattern formation but also in the fields of biology, 
mathematics and engineering. In the area of pool boiling so far only four studies of 
CML have been reported. They are Yanagita (1992), Shoji and Tajhna (1997), Shoji 
(1998) and Ghoshdastidar et al. (2004). So far no studies showing CML simulation of 
flow boiling have been found in the open literature. The present work is a first attempt 
in this direction. 

A CML is a dynamical system with discrete-time, discrete-space and continuous 
states. It usually consists of dynamical elements on a lattice which interact (are 
‘coupled’) with suitably chosen sets of other elements. The strategy of modelling 
dynamical phenomena in spatially extended systems by a CML model is based on the 
following steps. 

1 . Choose a (set of) macroscopic variables on a lattice. 

2. Decompose the processes imderlying the phenomena into independent components. 

3. Replace each'eomponent by a simple parallel dynamics on a lattice. 

4. Carry out each unit dynamics (or procedure) successively. 

In -the present study, the liquid is entering a constant wall temperature tube 
(Tw>Tsat) in a subcooled condition (Fig.l.l). The objectives of this study are to 
obtain qualitatively heat transfer coefficient, wall heat flux, mean fluid temperature, 
and vapour fraction as a function of axial coordinate. 

The modelling by CML is based on the assumption that the flow boiling is 
governed by (a) nucleation from cavities on the heated surface (b) forced convection 
and (c) phase chimge in the fluid bulk and mixing. The macroscopic variable chosen 



is temperature. The problem is transient and marches to steady state solution. The 
forced convection part in CML part is solved using CFD techniques, assuming 
average properties of the fluid. 

The reason why CML is qualitative is because the actual governing equations are 
not solved. Ideally, a computational fluid dynamics study of flow boiling would be 
possible if the Navier-Stokes equation could be solved together with continuity and 
energy equations along with phase change and appropriate boundary conditions. 
However, there are several problems in simulating the flow boiling phenomena based 
on CFD. First, numerical instability will occur because the boundary of the two 
phases is complex and varies with time. Second, the amount of CPU time needed for 
the integration of the equations will be very large since the calculations need to be 
carried out with high accuracy to avoid such numerical instabilities. Furthermore, 
with the numerical techniques available at present only a few bubbles can be dealt 
with at best. 

1.2 An overview of flow boiling 

Boiling in forced flow is called flow boiling or convective boiling. Generally, in most 
vapour-producing equipment vapourization occurs under forced convection. The flow 
conditions are influenced to a great extent by the pressure gradient along the heating 
surface. The vapour content increases along the path of flow up to the point of 
complete vaporization. In general, a liquid enters in a sub-cooled condition into a 
heated channel. Vapour bubbles formed at the wall due to nucleation condense again 
in the colder core of the liquid. If the liquid in the core is heated up to saturation 
temperature, then saturated nucleate boiling results. 

Figure 1.2 shows the various regimes in succession in the flow boiling in a 
vertical heated tube. After the bubbly flow characterizating nucleate boiling the 
individual bubbles grow together into large bubbles, that is, they develop into a slug 
flow. With increasing vapour content the large bubbles also grow together so that at 
first there is semi-annular flow and, subsequently, there forms at the tube wall a liquid 
film and a vapour core with liquid drops which is called dsmulsr-dispersed flow. With 
further addition of heat, the liquid film disappears downstream and in this part the 
two-phase flow contains vapour with liquid drops. This regime is called spray or mist 



flow. In most technical applications annular-dispersed flow occurs frequently. Slug 
flow occurs if the flow velocity is small. 

We already know that in pool boiling the flow field and heat transfer are 
determined by the difference between the hot surface and saturation temperature, as 
well as by the properties of the fluid and the heating surface. On the other hand, in 
boiling in forced flow the velocities of the vapoxir and liquid phases and the 
distribution of phases are additional factors that influence the flow field and heat 
transfer. Therefore, in flow boiling the empirical heat transfer coefficient relationships 
are of the form h =c (c{^) '' (m f(x), where m " and x* represent the mass flux of the 
fluid and the vapour quality respectively. In nucleate boiling the heat transfer 
coefficient is chiefly dependent on the heat flux and practically not at all on the flow 
velocity. On the contrary, in convective boiling the heat transfer coefficient is 
primarily influenced by the velocity of the flow or by the mass flux but hardly by the 
heat flux. In flow boiling, n is nearly zero, s lies between 0.6 and 0.8, while in the 
nucleate boiling portion, n = 0.75 and s lies between 0.1 and 0.3. 

For large Froude numbers, the flow regime is same both for vertical tube and 
horizontal tube, the reason being the forces of inertia are so great in comparison to the 
forces of gravity that no noticeable stratification occurs. With smaller Froude number, 
especially for Fr <.04, the stratification is very pronounced, and the influence of the 
slope of the tube must be considered. Fig. 1.3 shows various regimes of flow boiling 
in a heated horizontal tube for a non-stratified flow. 

There are two types of flow boiling crisis. One is film boiling and other is dryout. 
In film boiling crisis, the critical heat flux is reached at very high constant wall heat 
flux or very high constant wall temperature with very small vapotir quality. A vapom 
film forms at the wall that separates the liquid from the wall. Due to high thermal 
resistance of the vapour film, the wall heat flux drops and the critical heat flux is 
reached. In dryout boiling crisis, once the critical heat flux is reached, the liquid film 
disappears at the wall, which then becomes covered with vapour. In dryout the 
critical heat flux drops sharply with the quality. Dryout boiling occurs at high vapour 
quality. 



1.2 Literature Review 

Coupled map lattice (CML) has been recognized as a powerful tool of analysis to 
grasp the qualitative and fundamental nature of complex boiling phenomena and has 
been applied to many physical systems (Yanagita and Kaneko, 1993). The CML 
method is based on a dynamic system with continuous field variables but discrete 
space and time, in which local d 5 mamics propagates in space by diffusion or flow and 
time is advanced by repeated mapping. 

From the study of non-linear chaos dynamics it is known that a complex physical 
system is not always governed by a complex system of equations. Earlier studies of 
boiling based on non-linear dynamics, though few and recent, suggest that boiling is a 
kind of spatio-temporal chaotic phenomenon. The relevant papers in this regard are by 
Sadasivan et al. (1995), Shoji and Tajima (1997), Shoji (1998), Ellepola and Kenning 
(1996) and Nelson et al. (1996). However, applications to CML are not restricted to 
the problems in spatio-temporal chaos, but include pattern formations, some solid- 
state problems, biological information processing and engineering problems (Kaneko, 
1993). 

Using the CML method, Yanagita (1992) simulated the pool boiling 
phenomenon and succeeded in explaining the mode of transition from nucleate to film 
boiling. He assumed that fundamental dynamic processes of boiling are thermal 
convection, bubble rising motion and phase change. Although Yanagita’ s model is 
very attractive and path breaking, it deviates from actual boiling process at some 
points. Firstly, his model permits liquid to evaporate in the bulk, a phenomenon 
known as homogeneous nucleation, which is different from boiling in the usual sense. 
Secondly, in Yanagita’s model, pool boiling takes place even when the heater surface 
temperature is less than the saturation temperature of the liquid. This never happens in 
the actual boiling systems. Finally, the strict mode of film boiling is not realized by 
Yanagita’s model. 

Shoji (1998) corrected the deficiencies of the model of Yanagita by including 
nucleation sites on the heater surface and the Taylor instability. The model was 
applied to saturated and transient pool boiling of water on a small heated surface at 1 
atmospheric pressure. The effects of liquid subcooling and surface roughness were 
also investigated. Although Shoji’s model is basically sound, it has the following 
limitations. Firstly, the calculation of nucleation superheat distribution is sketchy and 



cannot be easily reproduced by other researchers. Secondly, no effect due to stirring 
action of bubbles is incorporated into his model. Finally, rather than using actual heat 
flux he used an apparent heat flux (without a unit) to plot the pool boiling curve and 
hence no information about the predicted CHF could be obtained. 

Ghoshdastidar et al. (2004) modified the basic theoretical model proposed by 
Shoji (1998) in terms of nucleation superheat distribution and mixing. The stirring 
action of the bubbles was modelled by increasing the fluid thermal diffusivity by an 
enhancement factor. The effectiveness of the enhancement factor approach in the 
model of Ghoshdastidar et al. (2004) is clearly seen in its capability of reproducing 
the saturated pool boiling curve well and predicting the critical heat flux (CHF) in the 
same order of magnitude of the actual value. 

So far no studies showing CML simulation of flow boiling have been found in the 
open literature. The present work is a first attempt in this direction. In this work non- 
stratified flow boiling of water as well as propane in a horizontal tube whose wall is 
maintained at constant temperature which is greater than the saturation temperature of 
the liquid at the entrance pressure is simulated by CML. The liquid is entering the 
tube in a sub-cooled condition. 


1.4 Objectives 

The objectives of present work are: 

(a) To develop a clear-cut methodology of producing a nucleation superheat 
distribution on the heated surface. 

(b) To find out the axial velocity distribution and mean velocity of the fluid at any 
circular cross section of the tube. 

(c) To find out temperature distribution and mean temperature of the fluid at any 
circular cross section of the tube. 

(d) To find out the mean heat transfer coefficient, mean heat flux and vapour 
fraction as a fimction of the axial coordinate. 

(e) To include the effect of stirring action of the bubbles. 

(f) To carry out a parametric study, for various levels of entry sub-cooling, wall 
temperatures and mass flow rates. 

(g) To find out the axial coordinate beyond which critical heat flux might occur. 
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Fig. 1.2 Various regimes 


of flow boiling in a vertical heated tube (Collier, 1972) 
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Fig 1.3 Various regimes of flow boiling in a heated horizontal tube for a non- 
stratified flow (Stephan, 1992) 


8 


CHAPTER 2 


PROBLEM FORMULATION 


2.1 Modeling 

In this chapter the problem formulation for the flow boiling of water is described but 
the general procedure is applicable to other fluids as well. Flow boiling of water at 
high pressure in a circular horizontal cylindrical tube is considered. The flow is 
assumed to be hydrodynamically fully developed but thermally undeveloped. 

In this work flow boiling is investigated in horizontal tube of diameter .02 m and 
length 9m. The temperature distribution in the tube wall is ignored as it has a very 
high thermal conductivity and small thickness. The three dimensional boiling field is 
approximated and modeling is based on CML. No prior assumption regarding the 
chaotic nature of the boiling phenomena needs to be made. 

2.2 Computational domain and lattices 

The computational domain in the r - 9 plane is divided into 20x20 lattices as shown 
in Fig. 2.1. Each lattice can contain either liquid or vapor. There are 20 grid points in 
the radial direction, 20 grid points in the angular direction and 10 grid points in the 
axial direction. The dots in Fig. 2.1 represent grid points. The broken lines in- the 
same Figure indicate the faces of lattices. Each lattice contains one grid point at its 
center. A typical lattice is shown by the shaded region in Fig. 2.1. Grid points at the 
boundary are surrounded by halfi-lattices, as shown in Fig. 2.1. A grid point is 
designated by (i, j, K) with ‘F increasing in the radial (r) direction, ‘y’ in the 
angular (^) direction and k in axial (z) direction. A lattice that encloses a grid point (z, 
y, k) is also referred to as {i,j, k). 



(i+lj, k) 



Fig. 2.1 Grid in (r, O') plane and pictorial representation of the treatment 
of the condition 
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2.3 Field variable 


For simplicity, temperature is employed as the only one field variable. In addition, a flag 
function, Fi, j, k is used for the convenience of calculations to show the phase of each 
lattice. Fi, j, k = ‘0’ and ‘T represent the lattice (i, j, k) in liquid and vapor phases 
respectively. 

2.4 Formulation of dynamic process 

In the CML method, dynamic processes are usually formulated in mappings. In the 
present model, it is assumed that boiling is governed by the following physics and 
dynamics. 

2.4.1 Nucleation on the heated surface 

According to Wang and Dhir (1993), nucleation cavities are distributed at random on the 
heated surface. Many cavities are distributed on each surface lattice but if it is assumed 
that every cavity has a conical shape, a larger cavity yields lower nucleation superheat. 
Therefore, only the cavity of maximum size [in this case, Dc as calculated firom equation 
(2.2)] is employed, since the active cavity of each surface lattice which determines the 
local nucleation superheat, ATact [equation (2.1)], is required for bubble nucleation. The 
term PR is added to Dm to create a randomness among the size of the large cavities. The 
nucleation superheat is given by 




Fig. 2.2 Circiunferential Nucleation Superheat Distribution for Water at the tube 


entrance 

Where Dc is the diameter of the largest nucleation cavity on a surface lattice. To 
calculate the following formula is used: 

Dc (m,y, k) =Dm+ |3R (2.2) 

Where Dc (m, j, k) represents the diameter of the largest nucleating cavity on the lattice 
(m,/, k), Dm is the minimum diameter of nucleating cavities on a heater surface lattice, (3, 
which is assigned a value of 0.99, indicates the maximum deviation from Dm and has 
same unit as Dc and Dm, i.e. micrometer, and R is a random number between 0 and 1 
assigned at the jth surface lattice. Although the number of grid points in the 
circtnnferential direction is taken as 20 in the present study, it may be varied. A total of 
1020 random numbers is generated using a NAG library random mnnber generator 
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subroutine. Every 50th random number of the set is employed at each surface lattice 
starting from the first. 

The minimum cavity diameter is arbitrarily taken as 1 |a,m. It is assumed that sites 
smaller than 1 pm do not exist on this particular heater surface. Thus, the maximum 
diameter of the nucleating cavities on a heater surface lattice becomes 1.99 pm as 
calculated from equation (2.2) when the random number takes on a value of 1 . 

Thus, using equations (2.1) and (2.2), the nucleation superheat distribution on the 
heated surface is calculated. Fig 2.2 shows a typical distribution of nucleation superheat 
at a particular axial location. 

2.4.1. 1 Nucleation Modelling 
Case I: (Nucleation not suppressed) 

If at a particular axial location (k) 

Tsat (k) “ Tmean (k) > 0, and if T (m, j. k) > Tact a k), where j = 1 to n 
and F (j, j, k) = 0, where i = ltom,j = lton 

then F (j, j, k) = 1 > where i = 1 8, j = 1 to n (2.3) 

if Tsat ■ Trnean^ 15 C. 

if Tsat - Tmean < 15 °C, then i = 10 in equation (2.3) 

Case 11: (Suppression of Nucleation) 

If at a particular axial location (k) 

Tsat (k) - Tmean (k) < 0, and if T (m, j. k) >Tact 0, k) and if F (j, j, k) = 0, where i = 1 to m, j = 1 to 
n, 

then F (m, j, k) = 1, where) = j3 to n (2.4) 

and F (i, j, k) = 1, where) = 1 to )3, i = 2 (2.5) 

To sum up the procedure of obtaining the nucleation superheat distribution, the following 
points should be kept in mind. The most important point in the calculation of nucleation 
superheat distribution is the minimum cavity diameter. Dm, on a surface lattice. It may be 
noted that the minimum cavity diameter is assumed to be the same for all surface lattices 
in the present problem. If, for a given surface. Dm is provided then the nucleation 
superheat distribution can be calculated by the combined application of equation (2.2) 
and equation (2.1). The following needs to be said about the physical (and mathematical) 
basis of equation (2.2). Since p is the maximum deviation from the minimum cavity 
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diameter, the product of (3 and the random number R, which takes a value between 0 and 
1 , essentially implies that a fraction of that maximum deviation is taken. The product 
[i.e. (3R ] is then added to the minimum value, Dm, to get Dc from equation (2.2). It can be 
easily seen that Dc can vary from 1 to 1.99 pm. The deviation P does not have to have 
fixed value. The choice of a value for P should be dictated by the information the analyst 
has about the minimum cavity size on a given surface lattice and the maximum cavity 
size on the heater surface. Thus, P may vary from one surface to other. 

2.4.2 Governing Equations: 

dv 

Assuming the flow to be laminar and fully developed (v,.= 0, — - = 0) and having no 

8z 

swirl (v g = 0), the z- momentum equation takes the form 


dv dp Id.dv 1 d V. 

p — ^ = — — + p (r — f- 

dt dz r dr dr r^ dO^ 


( 2 . 6 ) 


In the above equation, — is assumed to be a known quantity with p decreasing linearly 

dz 


withz. 

Att = 0, v_.= 0 
At r = 0, V . = finite 
At r == r 0 , V . = 0 (no slip) 

In addition the following periodic boundary conditions are valid. 


(2.7) 

( 2 . 8 ) 
(2.9) 


1 dv. 


r dO 
(V.-) 


r,e 


1 dv. 
^ 

r de\ 


r,$+27r 


r,0 


(V.-) 


r,d+27t 


The equation is; 

dT dT , 


dt 


dz 


11. (r^-) 

r dr dr 


+- 


kf d^T 


d9^ 


Att = 0,T = Tsat- ATsub 
At z = 0, T = Tsat- ATsub 
At r = 0, T = finite 
At r — r^, T — Tw (T^ ^ Tsat) 


( 2 . 10 ) 

( 2 . 11 ) 


( 2 . 12 ) 

(2.13) 

(2.14) 

(2.15) 

(2.16) 
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At = 10“^ s is used to cause stability. 

2.4.3 Accuracy: 

The solution has been validated with the analytical solutions for single phase flow and 
heat transfer (the fluid is water). Some representative plots are shown in Figure nos. 2.3, 

2.4 and 2.5. In Fig. 2.4 Nusselt number clearly approaches the steady state analytical 
solution, that is, 3.66 with increasing time. The analytical solution in Fig. 2.2 has been 
obtained from Bird et al. (1960). The analytical solution in Fig. 2.3 has been taken from 
Fox and McDonald (1995). 

2.4.4 Phase change and effects of bubble motion: 

It is known that latent heat is consumed when liquid evaporates and is released when 
vapor condenses. This phase change process is formulated as follows (Shoji, 1998): 


^ IJ,k ~ ^ f,J,k ^ ^c(i,j,k) 




. TP 


(2.19) 


If F ij,k = 1 -0 and T 


then ^ 


( 2 . 20 ) 


The suffix r|(i, j, k) represents the nearest neighboring four lattices and Tc(ij,k) is the phase 
change temperature, which is determined according to the phase change criteria given 
below. In equations (2.19) and (2.20), t) is a parameter related to the enthalpy of 
vaporization and has a unit of ‘^C. The aforesaid equations also represent the mixing 
effect of bulk liquid due to bubble motion. 

2.4.4. 1 Phase change criteria 

The value of Tc (i, j, k) for the liquid lattice adjacent to the wall is nucleation superheat, as 
shown in Fig.2.2. The value for liquid lattices in the bulk is assumed to be the 
homogeneous nucleation temperature of the liquid, but for liquid lattice neighboring the 
vapor lattice, is the saturation tenqjerature of the liquid. 
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Fig. 2.3 Unsteady state centre velocity in pipe flow 


Id 






Fig. 2.4 Fully developed steady laminar velocity profile in a pipe 




Fig 2.5 Transient Nusselt Number profile in a pipe flow 



CHAPTER 3 

METHOD OF SOLUTION 


3.1 Solution methodology 

In this chapter the method of solution by CML and overall solution algorithm is 
presented. 

3.1.1 Time advancement 

Time is advanced by repeating a set of mapping the dynamic processes (Section 2.4.2 
to 2.4.4) in such a manner that 

(3.1) 

Where the superscripts p and p+1 indicate the present and future time values of the 
temperature. 

3.1.2 Parameter values 

In the CML method, the parameter values are determined so as to reproduce the 
phenomena satisfactorily. In the present computations, the parameter value is 
ri = .01 

3.1.3 Overall solution algorithm 

1. Specify wall temperature Twaii = 230°C. 

2. Calculate fluid properties with pressure, as the pressure changes along the axial 
direction, due to pressure drop. 

3. Initialize temperature T, and velocity Vz- 

4. Calculate v ^ by solving momentum equation [equation (2.6)] . 

pi > 

5. Check Reynolds number. If Reynolds number exceeds 2000, then respecify — . 

oz 

6. Calculate by solving the energy equation [equation (2.12)]. 

7. Apply the nucleation criteria [equations (2.3), (2.4), (2.5)]. 

8. Determine the phase state of each fluid lattice by using phase change criteria. 

9. Apply equations (2.19) and (2.20) to obtain 

10. Check the phaise state of each fluid lattice again. 



11. Is steady state reached? If yes, go to step 13. If not calculate the overall vapour 
fraction, / (by volume), and equivalent density, specific heat and thermal conductivity 
and diffusivity, viscosity and kinematic viscosity of the hquid/vapour mixture using the 
following equations, enhance the molecular thermal conductivity and diffrisivity, 
viscosity and kinematic viscosity of the fluid by a suitable factor as described below to 
include the effect of mixing effect of the bubbles and go to step 4; 

keq =/kv+ (1-i) k/ 

Peq ~fPv'^ i^~J) Pi 

Ceq —fCy + (1-J) C I 

Meg=fMy+0--f)Ml 

The overall vapour fraction, /, is defined as the ratio of the volume of vapor and the 
volume of (liquid + vapour) mixture in the pool. 

For all values of f, the enhancement factor is 1.2. 

The enhancement factor basically takes into account the turbulent mixing in the fluid 

resulting from the bubble stirring action which enhances the thermal diffusivity, 
kinematic viscosity of the fluid. 

12. Calculate the steady state average nusselt number, average heat transfer coefficient a 
any circular cross section of the tube. 

13. From 13, calculate steady state wall heat flux. 
c(' (k) = hxnean (k)(Twaii - T^ean (k)), where k = 1 to 9. 

14. Print the axial location, heat flux, heat transfer coefficient, vapour fraction. 

15. STOP 

A concise flowchart of the solution algorithm is shown in Fig. 3.1 
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t= At 
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Fig 3. 1 Flow chart of the overall solution algorithm 











CHAPTER 4 


RESULTS AND DISCUSSION 


4.1 Introductory Remarks 

This chapter presents a parametric study of heat transfer coefficient, heat flux, fluid 
mean temperature and vapour fraction versus axial coordinate of the horizontal tube 
for various levels of entry subcooling, wall temperatures and mass flow rates. It may 
be noted that pressure changes along the axis and as a result the saturation 
temperature (Tsat) and thermophysical properties of the liquid and vapour also vary. 
Tsat and property variations with respect to pressure have been represented in the form 
of equations for convenience in computer implementation. These relations can be 
found in Appendix B for water as well as propane. The results shown here 
correspond to t = 1 sec. Note that A^ = 10'^s is used. Because of CPU time 
limitation, the results for higher time and steady state could not be obtained. The heat 
transfer coefficient actually is the mean value. It is defined as follows. 

2a- 


,(z) = 


A \q“(e,z)de 

In I _ q (z) 




(4.1) 


The heat flux is defined as 

ijt 


J! 


-j ZTT 


dT\ 


In J " dr 


\r=r^ 


de 


(4.2) 


The fluid mean temperature is basically the mixed mean temperature or mixing - cup 
temperature. 

1 r 1 


Av 


\vJdA 


2nro 


mean A 


^0 ^mean 0 0 


llvAr;9)nr,9)drdO 


(4.3) 


Where, 


. . 2!zr„ 

A j tiTq q q 


(4.4) 





The vapour fraction (f) is actually the volume fraction or void fraction which is defined 
as the ratio of the cross — sectional area (Ay) filled by the vapour at any location along the 
tube length to the total cross - sectional area (A) of the tube. Thus, 



These fractions do not change within a sufficiently small tube section Az . Thus, 


^ Az 
AAz 


(4.6) 


Therefore, the volume fraction (which is the same as the void fraction) of the vapour in 
the tube section under consideration is then 


f= 


V 


(4.7) 


4.2 Flow boiling of water 

4.2.1 hmeaii) <1 ) Tmeani f VS. Z for Twall ~ 230 C, ATjub” 30 C, Ptntrance~ 25 bar 

Figs 4.1 - 4.4 show hmean VS. z , vs. z, Tmean VS. z and f VS. z respectively for flow 

boiling of water when ^ = ~\.5Nlm\ T^aii = 230'’C and Tentrance= Tsa, at p - AT^^^, 

oz 

where A7’j„j= SO'^C. Fig. 4.1 shows that there is a steep drop of hmean near the entrance 

because of high Tw - Tsat- At z = 4m there is again a drop of heat transfer coefficient. 
This can be explained by the sudden rise of vapour fraction at that location as indicated in 
Fig. 4.4 and hence the conductivity of the liquid - vapour decreases because of higher 
vapour content. Thus the critical heat flux corresponding to z = 4m is around 2000 W/m^ 
as revealed in Fig. 4.2. Hence, z = 4 m may be considered as the length beyond which the 
tube runs the risk of getting irreversibly damaged. Fig. 4.3 indicates initial steep rise of 
mean temperature of the fluid followed by almost constant temperature region up to z = 4 
m. After this zone there is a drop in the mean temperature because of a drop in the heat 
transfer coefficient as revealed in Fig. 4.1. Tmean then further drops a little but shows a 
slight rise towards the end. This is because towards the exit of the tube vapour fraction is 
high and so is the vapour velocity resulting in more or less constant heat transfer 
coefiicient as the effect of low vapour conductivity is compensated by high vapour 
velocity. 
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4.2.2 Parametric Study 

4-2.2.1 hmean; ^ ? Tmean? f VS. z for different wall temperatures for ATsub ~ 

= -1.5 N/ml 

dz 

Fig. 4.5 shows that wall temperature has no effect on the heat transfer coefficient. 
Fig. 4.6 reveals that heat flux increases with wall temperature as Tw - Tsat is high. As 
expected greater wall temperature results in greater levels of mean fluid temperature 
as depicted in Fig. 4.7. Interestingly, increasing the wall temperature from 230°C to 
300°C does not change the vapour fraction (Fig. 4.8). This clearly shows that hmean is 
a function of vapour fraction. 

4.2.2.2 hmean; Q > Tmean; f VS. z for different ATsub’s for Twaii = 230®C and — = -1.5 

dz 

N/In^ 

Fig. 4.9 shows that with decrease in entry sub-cooling the heat transfer coefficient 
decreases in the section near the tube entrance but becomes independent of axial 
location in the latter part of the tube. Heat flux drops with lower levels of entry sub- 
cooling as revealed in Fig. 4.10, it being the lowest for 0°C, that is, when the 

fluid is entering at saturation condition. Fig. 4.1 1 shows that fluid mean temperature 
increases with decreased levels of entry sub-cooling. Fig. 4.11 also shows that as 
decreases the constant mean temperature zone increases. Fig. 4.12 depicts that 

near the tube entrance vapour fraction is high for low A Tsub’s whereas the maximum 
vapour fraction is same for all A Tsub’s- 

4.2.2.3 hmean; Tmean; f VS. z for different mass flow rates for Twaii = 230®C , 
ATsub = 30®C 

Fig. 4.13-4.16 describe the effect of fluid mass flow rate on variation of hmean, 
Tmean and f respectively along z. The mass flow rate of 0.348 x 10”^ Kg/s corresponds 

to — = -1.3 N/m^, while that of 0.401 xlO"^Kg/s corresponds to — = -1.5 N/m^ 
dz dz 

and that of 0.482 x 10'^ Kg/s corresponds to =-1.8 N/m^. The results in Fig. 4.13 

dz 

clearly show that heat transfer coefficient is virtually independent of the mass flow 
rate in the initial section of the tube where nucleate boiling effect is dominant. On the 
contrary, in the latter part of the tube convective vapourization is significant. Since 



with increased mass flow rates more vapourization occurs in the latter part of the tube 
(Fig. 4.16), heat transfer coefficient falls due to lower conductivity of vapour. Fig. 4.14 
vs. z) indicates the same trend as in Fig. 4.13. With increased mass flow the drop in 
the fluid mean temperature occurs at a lower axial location of the tube because of larger 
vapour formation (Fig. 4.15). Another interesting observation is that at the lowest mass 
flow rate of 0.348 x 1 Kg/s the rise of Tmean near the tube exit does not take place. As a 
matter of fact the mean temperature further falls. This is because the velocity of vapour 
is not high enough to compensate for the decrease in the vapour conductivity which 
results in lower heat transfer coefficient. 

4.3 Flow Boiling of Propane 

This section presents the results for flow boiling of propane. Fig. 4.17 shows a typical 
circumferential nucleation superheat distribution on the tube wall at the entrance. 

4.3.1 hinean» q y Tmcanj f VS. Z for Twall = 50®C, ATsub= 10®C, Pentrance = 5 bar 

Fig. 4.18 — 4.21 show hmean vs. z, vs. z, Tmean VS. z and f vs. z respectively for flow 

boiling of propane when-^ =.-3N/m\ Twaii = 50®C and Tentrance= Tsatatp,„,„„„ - 

OZ 

where AT^^^ = 30°C.The trends are similar to that for water except that the drop in heat 
transfer coefficient at some location down the tube is not as marked as that for water. 
Also, in Fig. 4.20 Tmean vs. z curve reveals that unlike in water boiling there is no 
constant temperature zone and there is a much larger latter section of the tube in which 
increase of temperature is visible. 

4.3.2 Parametric study 

4.3.2.1 1 hmean, q^ Tmean, f VS. z for different wall temperatures for ATsub = 10® C, ^ 

OZ 

= -.3N/m^ 

The relevant figures are Figs. 4.22 - 4.25. Here also similar trends are observed as for 
water. The heat transfer coefficient shows no dependence on wall temperature. 

4.3.2.2 hmean, q'^^, Tmean, f VS. z for different A Ts„b>s for Twaii = 40®C and ^=-.3N/m^. 

OZ 

The relevant figures are Figs. 4.26 - 4.29. The trends similar to water are observed. The 
exception is Fig. 4.28 where the fluid mean temperature is higher for ATsub = 
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20°C than for ATsub ~ 30°C in the section of the tube from z = 4.8 m to 6 m, the 
reason being higher vapour velocity in the former case due to larger vapour content in 
the convective vapourization zone which is larger than that for ATsub = 30'^C. 

4. 3. 2.3 hmean? 9 Tmean? f vs. z for different mass flow rates for Twaii = 40“C , A Ts^b 
= lO^C 

Figs. 4.30 - 4.33 describe the effect of fluid mass flow rate on variation of hmean, 
Tmean and f respectively along z. The mass flow rate of 0.292 x lO"^ Kg/s corresponds 

to — = -.1 N/m^, while that of 0.117x10"^ Kg/s corresponds to — = -.3 N/m^, and 
dz dz 

that of 0.175 X 10'’ Kg/s corresponds to — = -.6 N/m^. The results (Fig. 4.30) show 

dz 

that convective vapourization zone in each section is much larger as compared to the 
zone where nucleate boiling is dominant. As in the case of boiling of water, heat 
transfer coefficient is independent of mass flow rate while it decreases with increase 
in mass flow rate. The trend is same in Fig. 4.31 which depicts ((' vs. z variation. 
Fig. 4.32 shows that for m =0.292x10“^ Kg/s after a steep increase of mean 
temperature the gradient of temperature in the rest of the tube is constant but much 
lower than that near the entrance. It is also seen that generally Tmean increases for 
increase in mass flow rate except for one crossing of the temperature curves for in = 
0.117x10"’ Kg/s and tm =0.175x10"’ Kg/s. Fig. 4.33 reveals that vapour fraction 
remains constant throughout at a lowest mass flow rate, that is, w = 0.292x10"^ Kg/s 
showing very little effect of convective vapourization which for higher mass flow 
rates the effect of convective vapourization is clearly more visible. 
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Fig 4. 1 Heat Transfer Coefficient versus axial coordinate curve for Water 

^=-1.5N/m= 

dz 
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Fig. 4.3 Mean temperature of fluid (Tmean) versus axial coordinate curve for 

Water ^=-1.5 nW 
dz 
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Vapour Fraction 
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Fig. 4.5 Comparison of vs. z curves for different wall temperatures and 

at Ar . = 30°C , — =-1.5 N/m^ for Water 
dz 
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Fig.4.9 Comparison vs. z curves for different at constant wall 

temperature (Twall = 230 °C), ^ =-1.5 N/m^ for Water 
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Fig.4.1 1 Comparison of Tmean vs. z for different at constant wall temperature 

(Twall = 230 °C), — =-1.5 N/m^ for Water 
dz 



Vapour Fraction 




Heat T ransfer Coefficient (W/m2-K) 



Fig.4.13 Comparison of vs. z curves for different mass flow rates 

And at constant wall temperature (Twaii = 230 ° C), = 30*’C, 

for Water 
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Fig.4.14 Comparison of vs. z curves for different mass flow rates 
and at constant wall temperature (Twaii = 230 ® C), = 3 

for Water 
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Fig.4.15 Companson of Tn,ean vs. z for different mass flow rates 

and at constant wall temperature (T^„ = 230 °0 AT = 30“C 
for Water 
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Fig.4.16 Comparison of f vs z curve for different mass flow rates 
and at constant wall temperature (Twaii = 230 °C), = 


for Water 
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Grid point number in the circumferential direction 


Fig.4. 1 7 Circumferential Nucleation Superheat Distribution for Propane at the tube 
entrance 


43 



Heat Transfer C 


70 



44 




Heat Flux (W/m2) 



45 





Fig. 4.20 Mean temperature of fluid (Tmean) versus axial coordinate curve for 

Propane — =-.3 nW 
dz 
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Fig. 4.22 Comparison of vs. z curves for different wall temperatures and 

at = 10°C , — =-.3 nW for Propane 
dz 
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Fig. 4.24 Comparison of Tmean vs. z for different wall temperature at = 10®C, 

— =-.3 nW for Propane 
dz 
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Fig. 4.25 Comparison of f vs. z curves for different wall temperatures 
and at = 10°C, — =-.3 N/m^ for Propane 
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Fig.4.26 Comparison of vs. z cxirves for different at constant wall 

temperature (Twall = 40 “ C), — =-.3 N/m^ for Propane 

dz 
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Fig.4.27 Comparison of a[' vs. z curves for different at constant wall 

temperature (Twall = 40 °C), — =--.3 N/m^ for Propane 

dz 




frame 001 i 28 Dec 2004 j { | { 



Fig.4.28 Comparison of Tmean vs. z for different at constant wall temperature 

(Twall = 40 °C), =-.3 N/m^ for Propane 

dz 
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Fig.4.29 Comparison of f vs. z curve for different at constant wall temperature 

(Twall = 40 “ C), =-0.3 N/m^ for Propane 

dz 
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Fig.4.30 Comparison of vs. z curves for different mass flow rates 
and at constant wall temperature (Twaii = 40 ° C), = 1 0°C, 

for Propane 
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Fig.4.31 Comparison of vs. z curves for different mass flow rates 
and at constant wall temperature (Twaii 40 °C), =10 

for Propane 
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Fig.4.32 Comparison of Tmean vs. z for different mass flow rates 

and at constant wall temperature (Twaii= 40°C), 10°C, 

for Propane 
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CHAPTER 5 


CONCLUSIONS AND SCOPE FOR FUTURE WORK 

The present work shows the modelling of laminar, non - stratified flow boiling of 
water and propane in a constant wall temperature horizontal tube by the coupled map 
lattice method (CML). It may be noted that CML simulation produces only 
qualitative results. However, the simulation by CML is very useful to get an insight 
into complex phenomena such as boiling which can not be handled by available CFD 
techniques. The results are presented for variation of heat transfer coefficient, heat 
flux, mean fluid temperature and vapour fraction along the axial coordinate of the 
tube. The main observations from the results are the following, (i) The wall 
temperature does not affect the heat transfer coefficient and vapour fraction; (ii) The 
decrease in entry subcooling lowers the heat transfer coefficient near the tube entrance 
but does not influence the heat transfer coefficient in the latter part of the tube. 
Furthermore, vapor fraction is high near the tube entrance for low ’s whereas the 

maximum vapour fraction is same for all ’s. (iii) The heat transfer coefficient is 

virtually independent of mean flow rate in the initial section of the tube where the 
nucleate boiling effect is dominant. On the other hand, in the latter part of tlie tube 
where the convection vapourization is significant, because heat transfer coefficient 
falls with increase in mass flow rate. All these trends are similar for both water and 
propane. 

An exception for propane is noted where it is seen that fluid mean temperature is 
higher for AT^^^ = 20°C than for AT^^j = 30°C in the section of the tube from z = 4.8m 
to 6m. Another difference with respect to water is that the convective vapourization 
zone much larger as compared to the zone where the nucleate boiling is dominant. In 
addition to the above, the drop in heat transfer coefficient at some location down the 
tube is not as marked as that for water. 

The present analysis can be extended to include turbulent flow boiling (turbulent 
liquid, turbulent vapour) in a horizontal tube as turbulence is commonly encountered 
in industrial flow boiling. 
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APPENDIX A 


Finite difference based Method of Solution 
of Momentum and Energy Equations 


A.1 Handling of the condition at the centre 

Fig. A.l shows the grid in the r - 6 plane of the computational domain. Fig. A.2 
shows grid lines in the z — direction. 

Recall the z — momentum and energy equations, 
z - Momentum: 


dv dv dp \d.dv 1 


dt 

Energy: 

dt 


dz 


dz 


dT , 

pc„ h pc„v, — = kf 


r dr dr 


'1 d , dT' 

V — ) 

rdr dr 


de^ 


+ 


kf d^T 
r" d9^ 


(A.1) 


(A.2) 


It can be clearly seen from Eqs. (A.1) and (A.2) that 2'^'^ and 3’^^ terms of Eq. (A.1) and 
1®^ and 2"*^ terms of Eq. (A.2) become infinity at r = 0, leading to an undesirable 

situation. But, since at r = 0, neither nor — is necessarily zero, 1’ Hospital’s 

dr dr 


rule can not be applied. Furthermore, at r = 0, — - and — do not exist. Therefore, 

d0 de 

a special treatment at r = 0 is obviously needed. 

The aforesaid problem can be circumvented by taking a very small square region 
around the centre of the cylinder (see Fig. A.l) so that the governing differential 
equation in cartesian coordinates is valid in that region (Ghoshdastidar, 2004). The 
middle points of the top, bottom, right and left sides of the square will now coincide 
with the points on the horizontal and vertical lines passing through the centre of the 
circle (Fig. A.l). 

Therefore, at r = 0, i.e. at (1, 1, k), the z - momentum and energy equations take the 
following forms. 



z - Momentum; 


5v, 


dz 




d\,_ 

dx~ 


■ + 



Energy: 

dT 

PX^r. PCnV, 


dz 


a^r d^T 

dx^ dy^ 


(A.3) 


(A.4) 


Note Ax = Ay = A r in the square region. In the rest of the domain, x - 6 system is 
used. The method is quite accurate, as in a small region around the centre, a circle 
almost coincides with a square. Hence, the solution in the cartesian geometry is a 
good approximation. 


A.2 Discretization 

Equations (A.l) - (A.4) are discretized using Explicit finite - difference scheme 
(Ghoshdastidar, 2004). 

Discretizatized form of equation (A.l) 




- V 


zdJMi) 


At 






Az 




+ (a)o-./ 




.. - 2vf ; .. 


JP 




+ v: 


zii-\J,k) 




+ 




• V 






t+i) 


2Ar 
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{(^)(,vZi)f 


V 


P 

z(z,y+l,/:+l) 


2v 


p 

z(/, 7,A:+1) 


+ V 


P 

z{iJ-\,k-Pl) 


{aoY 


Valid for 


i = 2 m 

j = 1 n 

k = 1 1-1 


(A.5) 


Discretized form of equation (A.3) 


(p). 
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p+1 

^z(l,y,^+l) 




At 


+ (p)i 
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Az 
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r^y 




jP 

^2(2,1,A'+1) 


+ 

2(l,l,-t+l) ^ ^z(2J3,k+l) 


(Axf 

[;P — 

^2(2J2,A'+I) ^ ^z(2,j4Ml) 


+ - 




(A.6) 


Note that Ax = Ay = Ar. 
Discretized form of equation (A.2) 




(UMl) 


rpp+\ ^ rpp 
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rpp _ rpp 
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(/, 7,^+1) 
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1 
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Valid for 

i = 2 

j = l 

k= 1 


.m 

.n 

. 1-1 


Discretized form of equation (A.4) 


(S )(i, /„ 
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Note that Ax = Ay = Ar. 


(A.8) 
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A.3 Stability 

In order to ensure stability of the numerical solution At must be governed by the 
following equations. 


From equation (A. 5) 
At <■ 


V 


z{i,j,k+\) 


2(a)o-j, 






valid for i > 1 


(A.9) 


From equation (A. 6) 


At < 


"^^(1,1,^+!) 


2(a) 


p C 
(1, 1,^+1) 


^ (a)(u.*h-i) U^)' (Af)' 


■ +- 


(A.10) 


From equation (A. 7) 
At <T= 


V 





z(/,j,A'+l) 


, valid for i > 1 


(A.11) 


From equation (A.8) 


At< 


‘^2(1,1, i + l) 

tsz 


2(*/l 


f ^ 

f 2(1, 1,^+1) 


+ - 


P 2(l,l,/t+I) 


(Ax)^ (Ay) 


(A.12) 





X and • are almost coincident 



Fig. A. 1 Grid in (r, 9 ) plane and pictorial representation of the treatment 
of the condition 
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APPENDIX B 


Properties of Water and Propane as a function of pressure 


The following are the best fit relation for various thermophysical properties as a 
function of pressure for water and propane. The data have been taken from EES 
(Engineering Equation Solver) for Water and Perry’s Chemical Engineering 
Handbook for Propane. 

B.l Properties of liquid and vapour phases of water 

1. Surface tension: 

cr = (a + exp (-C (p — b)) 10”^ (B-1) 

Where, 

a = 2.4174, b = 0.0206, c = 192.5844 
cr is in N/m and p is in bar. 

2. Latent Heat of Vapoxuization: 

hj.^=a + bexp(-(-^^)^) (B.2) 

Where, 

a = -1216.0044, b = 2.933 x 10^ c = -4504.1362, d = 1731.561 
h is in KJ/Kg and p is in bar. 

3. Saturation Temperature: 

Tsat=ap^+bp^ + cp+ d (B.3) 

Where, 

a = 7.4677 x 1 0 b = -0.0302, c = 4.2744, d = 394.8 1 76 
Tsat is in Kelvin and p is in bar. 

4. Specific Heat of liquid and vapour: 
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Cp, = (a + exp(c(jr?-6))) (B.4) 

Where, 

a = 5.3514, b = 168.4687, c = 0.0938 
c ^ j is in KJ/Kg -K and p is in bar. 

c p,v = (a + exp (c {p-h)) (B.5) 

Where, 

a = 2.524, b = 52.9461, c = 0.0241 
c p ^ is in KJ/Kg -K and p is in bar. 

5 . Viscosity of liquid and vapour: 

,= (a + b exp ^ ) + e exp (-(^^— ^) " )10‘^ (B.6) 

d g 

Where, 

a = 65.2079, b = 1 143.1998, c = -15.409, d = 1 1.0574, e = 1 18.0526 
f= -81.0679, g= 127.1002 
; is in N-s/m^ and p is in bar. 


=(a + exp(c(p-b))) 

Where, 

a = 10.1242, b = -171.8871, c = 0.0084 

//„ is in N-s/m^ and.p is in bar. 

6. Density of liquid and vapour: 

1 . 

a -h b exp(-((p -c)ldf)\ O'" 


(B.7) 


(B.8) 


Wliere, 

a = 0.6558, b = 2.581 8, c = 483.0463, d = 351.4116 


68 



Pi is in Kg/m^ and p is in bar. 


Where, 

a = 0.0562, b - 1.8681, c = 0.4723 
is in Kg/m^ and p is in bar. 

7. Conductivity of liquid and vapour: 
k; =(a+bp)10”^ 

Where, 

a = 686.5352, b = -1.583 
k , is in W/m-K and p is in bar 

k^ = (a + bp)10'^ 

Where, 

a = 24.7355, b = 0.6088 
k^ is in W/m-K and p is in bar 


(B.9) 


(B.IO) 


(B.ll) 


In all the above expression the correlation coefficient varies between 0.9801 and 
0.998. 
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B.2 Properties of liquid and vapour phases of propane 

1 . Surface tension: 

cr = a+ b p + c exp (-p) 

Where, 

a = .01 13, b = -0.0004, = 0.1492 
<y is in N/m and p is in bar. 

2. Latent Heat of Vapourization: 

h^g=(a + bp + cp^) 

Where, 

a = 865.5829, b = 7.5152, c = -0.1596 
h is in KJ/Kg and p is in bar. 

3. Saturation Temperatiure: 

Tsat = (a + bp + cp^) 

(B.23) 

Where, 

a =234.3563, b = 8.1314, c = -0.1608 
Tsat is in Kelvin and p is in bar. 

4. Specific Heat ofliquid and vapour: 

Cp., = (a + bp + cp^) 

Where, 

a = 2.3097, b = 0.0364, c = 0.0007 
c p , is in KJ/Kg- K and p is in bar. 

=(a + bp) 


(B.21) 


(B.22) 


(B.24) 


(B.25) 
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Where, 


a = 1.5294, b = 0.0613 
c p ^ is in KJ/Kg-K and p is in bar. 

5. Viscosity of liquid and vapour: 

// ,= (a + b exp (-p) + cp) 10"^ (B.26) 

Where, 

a= 1.5172, b = 1.9121, c = -0.04 
, is in N-s/m^ and p is in bar. 

//„=(a + bp + cp^) 10*^ (B.27) 

Where, 

a = 0.0730, b = 0.0012, c = 0.000013 

is in N-s/m^ and p is in bar. 

6. Density ofliquid and vapour: 

Pj = (a + b p + c exp (-p)) (B.28) 

Where, 

a = 553.8787, b = -6.0358, c = 122.95 
Pi is in Kg/m^ and p is in bar. 

p„(k) = (a + bp) (B.29) 

Where, 

a = -2.0782, b = 2.4746 
p^ is in Kg/m^ and p is in bar. 
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6. Conductivity of liquid and vapour; 
k,= (a + b p + c exp (-p)) 


(B.30) 


Where, 

a = 0;i098, b = -0.0016, c = 0.0762 
k, is in W/m-K and p is in bar 

k„ =(a-Hbp + cp^ ) (B.31) 

Where, 

a = 0.01 15, b = 0.00101, c = -0.000015 
k„ is in W/m-K and p is in bar 

In the all above expressions the correlation coefficient various between 0.98 and 
0.998. 
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